{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 移动平均(MA)模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 原理讲解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 模型介绍\n",
    "\n",
    "另一类时间序列模型为“移动平均过程”（Moving Average Process，简记 MA）。记一阶移动平均过程为 MA（1）：\n",
    "\n",
    "$$y_{t}=\\mu+\\varepsilon_{t}+\\theta \\varepsilon_{t-1}\\tag{1}$$\n",
    "\n",
    "其中，{$\\varepsilon_{t}$} 为白噪声，而 $\\varepsilon_{t}$ 的系数被标准化为1。由于 $y_t$ 可以被看成是白噪声的移动平均，故得名。考虑使用条件 MLE 估计。为此，假设 {$\\varepsilon_{t}$} 为独立同分布且服从正态分布 $N\\left(0, \\sigma_{\\varepsilon}^{2}\\right)$。如果已经知道 {$\\varepsilon_{t-1}$} ，则\n",
    "\n",
    "$$y_{t} | \\varepsilon_{t-1} \\sim N\\left(\\mu+\\theta \\varepsilon_{t-1}, \\sigma_{\\varepsilon}^{2}\\right)\\tag{2}$$\n",
    "\n",
    "假设 $\\varepsilon_0=0$，则可以知道 $\\varepsilon_{1}=y_{1}-\\mu$。给定 $\\varepsilon_{1}$，则可以知道 $\\varepsilon_{2}=y_{2}-\\mu-\\theta \\varepsilon_{1}$。以此类推，则可以使用递推公式“$\\varepsilon_{t}=y_{t}-\\mu-\\theta \\varepsilon_{t-1}$”计算出全部 $\\left\\{\\varepsilon_{1}, \\varepsilon_{2}, \\cdots, \\varepsilon_{T}\\right\\}$。这样，在给定 $\\varepsilon_0=0$ 的条件下，就可以写下样本数据 $\\left\\{y_{1}, y_{2}, \\cdots, y_{T}\\right\\}$ 的似然函数，然后使用数值方法求解其最大化问题。\n",
    "\n",
    "更一般地，对于 $q$ 阶移动平均过程，记为 $MA（q）$：\n",
    "\n",
    "$y_{t}=\\mu+\\varepsilon_{t}+\\theta_{1} \\varepsilon_{t-1}+\\theta_{2} \\varepsilon_{t-2}+\\cdots+\\theta_{q} \\varepsilon_{t-q}\\tag{3}$\n",
    "\n",
    "也可以进行条件 MLE 估计，即在给定“$\\varepsilon_{0}=\\varepsilon_{-1}=\\cdots=\\varepsilon_{-q+1}=0$”的条件下，最大化样本数据的似然函数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MA(1) 模型最大似然函数推导\n",
    "\n",
    "$y_t=μ+ε_t+θε_(t-1)$    \n",
    "\n",
    "$y_1=μ+ε_1+θε_0=μ+ε_1  ⇒ε_1=y_1-μ  ⇒y_1∼N(μ,σ_ε^2)$   \n",
    "\n",
    "$y_2=μ+ε_2+θε_1=(1-θ)μ+θy_1+ε_2  ⇒ε_2= y_2-θy_1-(1-θ)μ  ⇒y_2∼N((1-θ)μ+θy_1,σ_ε^2)$     \n",
    "\n",
    "以此类推：$y_3∼N((1-θ+θ^2)μ+θy_2-θ^2 y_1,σ_ε^2)$       \n",
    "\n",
    "因此：\n",
    "\n",
    "$$f\\left(y_{1}\\right)=\\frac{1}{\\sqrt{2 \\pi \\sigma_{\\varepsilon}^{2}}} e^{-\\frac{\\left(y_{1}-\\mu\\right)^{2}}{2 \\sigma_{\\varepsilon}^{2}}}, f_{y_{2} | y_{1}}\\left(y_{2} | y_{1}\\right)=\\frac{1}{\\sqrt{2 \\pi \\sigma_{\\varepsilon}^{2}}} e^{-\\frac{\\left(y_{2}-(1-\\theta) \\mu-\\theta_{1}\\right)^{2}}{2 \\sigma_{\\varepsilon}^{2}}}\\tag{4}$$    \n",
    "\n",
    "$$\\ln L=-\\frac{1}{2} \\ln 2 \\pi-\\frac{1}{2} \\ln \\sigma_{\\varepsilon}^{2}-\\frac{\\left(y_{1}-\\mu\\right)^{2}}{2 \\sigma_{\\varepsilon}^{2}}-\\frac{T-1}{2} \\ln 2 \\pi-\\frac{T-1}{2} \\ln \\sigma_{\\varepsilon}^{2}-\\sum_{t=2}^{T} \\frac{\\left(y_{t}-c_{1} \\mu-c_{2}\\right)}{2 \\sigma_{\\varepsilon}^{2}}\\tag{5}$$   \n",
    "\n",
    "其中：$c_1=1-θ+θ^2-θ^3+......;    c_2=θy_(t-1)-θ^2 y_(t-2)+θ^3 y_(t-3)+......$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## statsmodels 库实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matlab实现\n",
    "\n",
    "MAlnL.m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```matlab\n",
    "function lnL = MAlnL(Beta,Y)\n",
    "%Beta(1) = u ;Beta(2) = sita;Beta(3) = sigma2;\n",
    "T = length(Y);\n",
    "lnL1 = -0.5*log(2*pi)-0.5*log(Beta(3)^2)-(Y(1)-Beta(1))^2/(2*Beta(3)^2)-0.5*log(2*pi)*(T-1)-0.5*(T-1)*log(Beta(3)^2);\n",
    "for i = 2:T\n",
    "    for j = 1:i\n",
    "        c11(j) = (-Beta(2))^(j-1);\n",
    "    end\n",
    "    c1 = sum(c11)*Beta(1);\n",
    "    for j = 1:i-1\n",
    "        c22(j) = Y(i-j)*(Beta(2)^j)*(-1)^(j+1);\n",
    "    end\n",
    "    c2 = sum(c22);\n",
    "    lnL2(i) = (Y(i)- c1 - c2)^2;\n",
    "end\n",
    "lnL = lnL1-0.5*sum(lnL2)/Beta(3)^2;\n",
    "lnL = -lnL;\n",
    "\n",
    "% %导入数据\n",
    "% data = xlsread('C:\\Users\\Administrator\\Desktop\\hourse.xlsx');\n",
    "% f1 = data(:,2); f2 = data(:,3); e = data(:,6);\n",
    "% \n",
    "% % 2.初始参数设定\n",
    "% maxsize         = 2000;         % 生成均匀随机数的个数(用于赋初值)\n",
    "% REP\t\t\t    = 100;          % 若发散则继续进行迭代的次数\n",
    "% nInitialVectors    = [maxsize, 3];    % 生成随机数向量\n",
    "% MaxFunEvals    = 5000;         % 函数评价的最大次数\n",
    "% MaxIter         = 5000;         % 允许迭代的最大次数\n",
    "% options = optimset('LargeScale', 'off','HessUpdate', 'dfp','MaxFunEvals', ...\n",
    "% MaxFunEvals, 'display', 'on', 'MaxIter', MaxIter, 'TolFun', 1e-6, 'TolX', 1e-6,'TolCon',10^-12);\n",
    "% \n",
    "% % 3.寻找最优初值\n",
    "% initialTargetVectors = unifrnd(0,10, nInitialVectors);\n",
    "% RQfval = zeros(nInitialVectors(1), 1);\n",
    "% for i = 1:nInitialVectors(1)\n",
    "%     RQfval(i) = MAlnL (initialTargetVectors(i,:), f1);\n",
    "% end\n",
    "% Results          = [RQfval, initialTargetVectors];\n",
    "% SortedResults    = sortrows(Results,1);\n",
    "% BestInitialCond  = SortedResults(1,2: size(Results,2));    \n",
    "% \n",
    "% % 4.迭代求出最优估计值Beta\n",
    "% [Beta, fval exitflag] = fminsearch(' MAlnL ', BestInitialCond,options,f1);\n",
    "% for it = 1:REP\n",
    "% if exitflag == 1, break, end\n",
    "% [Beta, fval exitflag] = fminsearch(' MAlnL ', BestInitialCond,options,f1);\n",
    "% end\n",
    "% if exitflag~=1, warning('警告：迭代并没有完成'), end\n",
    "\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
